Revealing Superfluid— Mott-Insulator Transition in an Optical Lattice 
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We study (by an exact numerical scheme) the single-particle density matrix of ~ 10 s ultracold 
atoms in an optical lattice with a parabolic confining potential. Our simulation is directly relevant 
to the interpretation and further development of the recent pioneering experiment jl| . In particular, 
we show that restructuring of the spatial distribution of the superfluid component when a domain 
of Mott-insulator phase appears in the system, results in a fine structure of the particle momen- 
tum distribution. This feature may be used to locate the point of the superfluid-Mott-insulator 
transition. 



The fascinating physics of the superfluid-insulator 
transition in a system of interacting bosons on a lattice 
has been attracting constant interest of theorists during 
recent years [||-^). Lattice bosons is one of the simplest 
many-body problems with strong competition between 
the potential and kinetic energy and a typical exam- 
ple of the quantum phase transition system. One of its 
great advantages is the possibility to study it by powerful 
Monte Carlo methods which nowadays allow simulations 
of many thousands of particles at low temperature with 
unprecedented accuracy (see, e.g., ||). However, until 
very recently the canonical Bose-Hubbard model 
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(where a\ creates a particle on the site i, <ij> stands for 
the nearest-neighbor sites, ni = ajoj, and t, U, and /i,, 
are the hopping amplitude, the on-site interaction, and 
the on-site external field, respectively) was not particu- 
larly useful in the analysis of realistic systems. The situ- 
ation has changed with the exciting success of the exper- 
iment by Greiner et al. (originally proposed by Jaksch 
et al. ]lC[| ) in which a gas of ultracold 87 Rb atoms was 
trapped in a three-dimensional, simple-cubic optical lat- 
tice potential. The uniqueness of the new system is that 
it is adequately described by the Bose-Hubbard Hamil- 
tonian ifToUTl, and allows virtually an unlimited control 
over the strength of the effective interparticle interaction 
U/t and particle density. 

The characteristic feature of the experimental setup of 
Rcf. jl| is the presence of the overall parabolic potential 
V(r) which confines the sample. This feature could be 
of great advantage if one was able to directly measure 
the spatial density distribution in the trap. We recall 
the structure of the fi — U/t phase diagram for the Bose- 
Hubbard system [Q which predicts commensurate parti- 
cle density distribution for the insulating phase whenever 
the chemical potential lies within the Mott- Hubbard gap. 
The slowly varying (at the length scale of the lattice pe- 



riod) trapping potential effectively provides a scan over 
fj, of this phase diagram at a fixed value of U/t. 

Unfortunately, what is measured in the experiment is 
not the original spatial density distribution in the trap, 
but the absorption image of the free evolving atomic 
cloud, after the trapping/optical potential is removed; i.e. 
the quantity which is directly related to the the single- 
particle density matrix in momentum space, pkk = ^k- 
[This statement implies that in the free evolving atomic 
cloud the interparticle interaction can be neglected; see 
the discussion below.] Now, in terms of rik the inhomoge- 
neous trapping potential is a disadvantage since it broad- 
ens the superfluid 8- functional contribution at k = 0, 
and the observed picture is a convolution of the original 
real-space density matrix /?(r, r'). As we show below one 
has then to look at the fine-structure of the central peak 
in the experimental data to decipher the Mott-Hubbard 
phase diagram (the "fading" of the Bragg peaks in the ex- 
periment has very little to do with the superfluid-Mott- 
insulator transitions and happens when the system is al- 
ready deep in the insulating phase). 

In this Letter we relate quantitatively the particle dis- 
tribution in momentum space observed in experiments 
to the corresponding spatial density distribution in the 
trap. Our ultimate goal is to reveal which features (if 
any) in the structure of indicate unambiguously the 
presence of the Mott phase. To this end we perform 
quantum Monte Carlo simulations of the single-particle 
density matrix for the Bose-Hubbard systems with up to 
16 3 lattice sites using continuous-time Worm algorithm 
ifTTf . We find that the onset of the phase transition in 
the trap center should result in appearance of at least 
one satellite peak in nk, reflecting a shell- type form of 
the superfluid component. This peak was not mentioned 
in the experiments of Ref. |Q. We suggest a possible 
explanation to this fact, and argue that by collimating 
the expanding atomic cloud one can render this peak ob- 
servable. We also discuss the role of self-repulsion in the 
free expanding cloud, which can affect the simple inter- 
pretation of the absorption images in terms of the initial 
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single-particle momentum distribution. 
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FIG. 1. Particle density distributions (on-site filling factors) as functions of the lattice site distance from the trap center for 
various coupling parameters and filling factors in the center: U/t = 24 (a), U/t = 32 (b), U/t = 80 (c), U/t — 80 (d), U/t — 80 
(e), E7/t = 80 (f). 



In Figs, l(a-f) we present our data for the density dis- 
tribution as function of the lattice cite distance from the 
trap center, r/a, where a is the lattice constant. For all 
practical purposes one may assume the zero-temperature 
limit here. The simulation was done at finite but very low 
T = 0.2t; the relevant energy parameters in this model 
arc the bandwidth W = I2t and (U/t) c ~ 35 In 
accordance with the phase diagram of Ref. [Q], we ob- 
serve a shell-type structure of the particle density with 
the Mott-insulator phases visible as integer plateau re- 
gions. 

Next, we relate each of the above figures to the cor- 
responding momentum distribution function n^. By 
definition, — J d 3 rd 3 r' exp[ik(r — r')]p(r, r'), where 
p(r, r') = (V^(r)i/>(r')}, and ^(r) is the bosonic field op- 
erator. In our case of a single-zone lattice, the field op- 
erator is expanded as follows: 

<Kr) =5>(r- ri K , (2) 

i 

where <j> is the Wannier function. We thus finally have 



n k =\ < P(k)\ 2 J^e ik ^-^p ij , (3) 

where 

Pij = ( a i a j) . ( 4 ) 

and </>(k) is the Fourier transform of 4>(r). From Eq. (0) 
it is seen that up to a trivial reweighting factor |</>(k)p 
the distribution is a periodic function in the reciprocal 
lattice. Thus without loss of generality we may restrict 
ourselves to the first Brillouin zone. Actually, (/>(k) has 
nothing to do with the Bose-Hubbard model, being a 
non-universal property of the lattice cite potential; in 
what follows we will ignore this function altogether by 
formally setting it to unity. 

Having calculated pij with the Worm algorithm 
we readily obtain using Eq. (^|); the results are pre- 
sented in Figs. 2(a-f). In Fig. 2(a) we see a typical pic- 
ture for the strongly correlated superfluid phase, char- 
acterized by a single, narrow peak at small momenta. 
When a domain of the Mott-insulating phase appears 
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in the center of the trap (where the on-site filling is 
close to unity), a pronounced fine structure develops in 
Fig. 2(b). We associate this structure with the shell- 
type form of the condensate wave-function. To prove the 
point we model the situation with the pure-condensate 
density matrix p(r,r') = ^(r)^^), where the conden- 
sate wave- function 5'o( r ) has the shell- type form with 
the shell radius I. The presence of the Mott insulator is 
taken into account through the suppression of the ^o(r) 
in the center. The Fourier transform of such 5'o( r ) is 
alternating in sign, with the half-period in k related to 



the shell radius as k ~ ir/l. Thus in the pure condensate 
we would see exact zeros in nk with the typical separa- 
tion between them ~ ir/l. Surprisingly, this naive model 
works extremely well and adequately describes the case 
of the realistic strongly correlated system close to the 
phase transition (in Fig. 2(b) the coupling U/t — 32 is 
close to the critical value estimated in Refs. (||^]). We 
consider the appearance of the satellite peaks as a clear 
signature of the Mott-insulator transition in the center 
of the trap. 





FIG. 2. The nk distributions in the first Brillouin zone in 
matrices for systems shown in Figs, l(a-f). 

In Fig. 2(c) the coupling strength is significantly in- 
creased to U/t = 80, but the shell- type structure of the 
superfluid phase, and the corresponding fine-structure 
in rT-k, is still present. The crucial difference between 
Fig. 2(c) and Fig. 2(d) is in the suppression (almost com- 
plete) of the superfluid fraction. Now the distribution rik 
has only the central peak and is extended towards large 
momenta, as expected for the system deep in the insu- 
lating phase. Still, it is not flat (!), which tells us about 
large off-diagonal correlations between the nearest lattice 
cites even for U/t as large as ~ 2.5 (U/t) c . 

Fig. 2(e) is similar in physics to the case (a), except 



the (0, 0, 1) direction derived from the single-particle density 

for the large-momentum tail due to the Mott-insulator 
shell. Finally, in Fig. 2(f) we again see the fine structure 
of satellite peaks reflecting the appearance of the Mott- 
insulator phase in the center of the trap and the cor- 
responding superfluid phase shell [in close resemblance 
with Fig. 2(b)]. 

We note, that momentum distributions fik presented 
above, may be observed experimentally if atoms are col- 
limated out of the expanding cloud so that the distribu- 
tion in a given direction is photographed. In the cur- 
rent setup, the absorption images of the three dimen- 
sional distribution are taken along two orthogonal axis. 
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This procedure reveals only the integral N(k x ,k y ) cx 
J^°oo dk z n{k). It is easy to see that integration effectively 
erases fine-structure features of n(k) — although peaks do 
not disappear completely, they now show up as shoulders 
in N(k x , k y ). 

Finally, we would like to discuss the problem of the 
repulsive interaction between particles during the initial 
period of their free expansion. Obviously, the interpreta- 
tion of the final absorption pattern in terms of the initial 
momentum distribution is valid only if the effect of this 
interaction is small enough. Meanwhile, given realistic 
experimental parameters, this turns out to be the case 
only for rather moderate system sizes. The criterion for 
neglecting the effect of interparticle repulsion is 



-Ekin/^pot ^> 1 , 



(5) 



where -Ekin and E pot are, respectively, the kinetic and 
potential energy per particle in the most fragile low- 
momentum part of the distribution rik and at the most 
dangerous period of free evolution at the end of the 
restructuring period, when the "discrete" distribution 
of density transforms into the "continuous" spatial dis- 
tribution with the typical size of order of the original 
system size (plus the corresponding replicas in higher 
Brillouin zones). For the potential energy we have 
Epot ~ nllo, where n is the continuous number density, 
Uq = ATrh 2 a s /m, a s is the s-scattering length, and m is 
the atom mass. Recalling that the lattice filling factor is 
of order unity, we can estimate n ~ 1/a 3 . The lowest ki- 
netic energy is associated with the spatial distribution of 
the condensate. Estimating ~ it 2 Ti 2 /ma 2 L 2 , where 
integer L stands for the typical size of the superfluid com- 
ponent in units of lattice constants, we arrive at a simple 
requirement 



(a/a s ) L 2 > 1 . 



(6) 



In the experiment of Ref. Q], the ratio a/a s is of order 
10 2 . Hence, we are restricted to L < 10, that is to typical 
system sizes of our present simulation. Note also, that 
the condition (||) is much easier to satisfy for the atomic 
cloud in the second Brillouin-zone peak, where the spatial 
density is significantly suppressed [Q . 

Summarizing, we presented a simulation of the ground 
state properties of ultracold atoms in an optical lattice 
with confining external potential, in the regime where 
the Mott-insulator and superfluid phases coexist. We 



have demonstrated that when the insulator domain in the 
center of the trap is surrounded by the superfluid com- 
ponent, the global momentum distribution of particles 
features satellite peaks. This picture can be employed by 
the experiment as an unambiguous evidence of the Mott 
transition. We do not see other features of the momen- 
tum distribution that could be associated with the Mott 
transition: unless the ratio U jt is not much larger than 
the critical one, the momentum distribution in the recip- 
rocal lattice still has a peaked form reflecting strong local 
off-diagonal correlations. 
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